Add k-infinity estimator for IndependentOperator - #4053

Looks like a nice PR.
I have some structural suggestions:

1. Calculate kinf by default when there is some fission reaction and every fission score has a corresponding nu-fission score.
   That way we do not need the calculate_kinf keyword argument.
2. Extend chain.REACTIONS to save all emitted particles. Then use it to figure out how many neutrons are emitted for every supported reaction instead of using regex.
3. Extending Nuclide::collapse_rate to collapse nu_total according to the fission * flux spectrum. That way we plug into the existing machinery instead of adding a special python function to handle the nu-fission case.

Thanks for the thorough review @GuySten! I've addressed all three suggestions:

1. Auto-detection instead of calculate_kinf parameter -- Removed calculate_kinf from __init__ and from_nuclides. The operator now automatically estimates k-infinity when keff is not provided and all MicroXS instances contain both fission and nu-fission reactions.

2. neutrons_out in REACTIONS -- Extended ReactionInfo with a neutrons_out field for all 84 reactions in chain.py. This replaces the _XN_REACTION regex and _neutrons_emitted() function. The loss calculation in _estimate_k_inf now uses REACTIONS[rxn].neutrons_out directly.

3. nu-fission collapsing in C++ -- Moved the nu-fission cross section collapsing from Python to Nuclide::collapse_nu_fission_rate in C++, with a C API and ctypes binding. This follows the same temperature interpolation and trapezoidal integration pattern as collapse_rate, but integrates v(E)*sigma_f(E)*flux(E).

All existing kinf tests updated and passing.

Thanks for the rework. The neutrons_out table is a real improvement over the
regex, the C++ move is the right call, and the auto-detection is in the right
place. I've verified the table against OpenMC's own convention (§4) and it's
exactly right.

I've pushed a commit with the mechanical cleanups (strict REACTION_MT lookup,
version markers, an empty-micros guard, a warn() on zero loss, N_FISSION,
and a clang-format fix that was failing the C++ Format Check). The rest is
below. The first three items are about what the estimator computes, and they're
upstream of the two code-level ones.

1. What quantity is this?

The docstring says "the infinite multiplication factor," which is a statement
about the system. What the code computes is a neutron balance over
self.local_mats, the depletable materials only. Non-depletable moderator,
clad, structure and reflector contribute exactly zero to the loss term, and
there's no data path by which they could: _get_burnable_mats keeps only
mat.depletable (openmc_operator.py:213), so the operator holds no number
densities and no MicroXS for anything else.

For a single homogenized material, or a model where everything is marked
depletable, those coincide. For a pin cell with non-depleting water and Zr clad
they don't. The one-group absorption fraction in fuel is 1 by construction, so
the estimate is high by roughly 1/f where f is that fraction. For a fresh UO2
pin cell that is a few percent, and considerably more with soluble boron in the
moderator, drifting with burnup as the fission-product and boron balance
changes. Because the estimate is enabled automatically and
silently, a user in that situation gets a plausible-looking number out of
Results.get_keff() with nothing marking it approximate.

I'd like the docstring to state the estimator explicitly, with the index ranges,
since that's where all the content is:

$$k_\infty = \frac{\sum_m \frac{1}{V_m}\sum_i N_{m,i}\sum_g (\nu\sigma_f)_{i,g},\phi_{m,g}} {\sum_m \frac{1}{V_m}\sum_i N_{m,i}\sum_g\sum_j (1-x_j),\sigma_{i,j,g},\phi_{m,g} ;+; \Lambda_\mathrm{bg}}$$

• m: depletable materials only
• i: nuclides in both the supplied MicroXS and the number-density array
• j: reactions present in the supplied MicroXS
• φ: the user-supplied flux, frozen across the depletion
• x_j: neutrons emitted by reaction j
• Λ_bg: absorption outside the depletable materials, §3

Worth stating plainly in the docstring that this is k∞ and not k-eff. The two
are related by k∞ = k-eff / (1 - L), so k∞ exceeds k-eff whenever there is
leakage and they coincide only when L = 0.

2. A narrowed reactions list silently reports nu-bar

This one is independent of §1 and §3, and I think it should be fixed regardless
of what you decide there.

'nu-fission' is never in chain.reactions, which is built from the chain
file's transmutation reaction types. So the auto-detection can only ever fire
for a user who passed an explicit reactions list. That means the at-risk
population is exactly the population hand-writing that list.

Nothing currently stops reactions=['fission', 'nu-fission']. That satisfies
the auto-detect, 'nu-fission' goes to production, 'fission' falls through to
loss because there's no 'fission' key in REACTIONS, and nothing else
contributes:

k = νΣf / Σf = ν̄ ≈ 2.43

reported silently as the eigenvalue. That exact list is a bit contrived, but the
realistic version, dropping one or two absorption channels from an otherwise
sensible list, fails the same way with a less obvious number.

The cheap fix is a static precondition at construction rather than a tally:

set(self.chain.reactions) <= set(m.reactions) and 'nu-fission' in m.reactions

self.chain.reactions is the right reference because it's the same list
_calculate_reaction_rates indexes the rate array with
(openmc_operator.py:477). If the MicroXS covers every reaction the operator
will deplete with, it covers every absorption channel the chain knows about.

Two plumbing notes. self.chain doesn't exist yet where _calculate_kinf
currently sits (independent_operator.py:149), since super().__init__() is
the last statement and the chain loads in abc.py:154, so the check has to move
below the super call. That's also correct on the merits, because
reduce_chain_level rewrites self.chain at openmc_operator.py:138 and you
want to check against the reduced chain.

A companion check on nuclides is worth having for the same reason: every chain
nuclide with cross-section data should appear in m.nuclides, since
from_multigroup_flux silently skips nuclides without data.

This can't see gaps in the chain itself, a reduced chain missing Xe-135 and
Sm-149 for instance. That's a "your chain is inadequate for this problem" issue
rather than a misconfiguration, and I'd cover it with a sentence in the
docstring rather than machinery.

3. Λ_bg, the absorption outside the depletable materials

Every neutron either is absorbed or leaks, so per source neutron in an
eigenvalue run:

$$A + L = 1 + X$$

with A absorption (disappearance + fission), L leakage, X net (n,xn) production
Σ(x-1)R. For k∞ we want P/(A - X), and the identity gives A - X = 1 - L
directly. X cancels, so it never needs tallying, in the depletable region or
in the background. Splitting into the part the operator recomputes each step and
the part it can't:

$$(A_\mathrm{dep} - X_\mathrm{dep}) + \Lambda_\mathrm{bg} = 1 - L$$

The first bracket is already your loss, since (1 - x_j) gives you A - X in
one pass. So:

1. Read L and k-eff off the statepoint (StatePoint.global_tallies) and return
   them from get_microxs_and_flux.
2. At BOL compute loss_BOL from the MicroXS and number densities.
3. Λ_bg = (1 - L) - loss_BOL, one scalar, stored once.
4. Each step: k_inf = production(t) / (loss(t) + Λ_bg).

No new tallies, no new scores, no transport changes. The only API change is two
floats coming back from a function that already runs the solve. Worth doing that
compatibly rather than widening the return tuple, since get_microxs_and_flux
is public and every existing caller unpacks two values. Attaching them to the
returned MicroXS objects, or an opt-in keyword, both avoid the break.

Background materials don't change composition, so freezing Λ_bg is the same
assumption the operator already makes about the flux and the cross sections, and
introduces no new approximation. An all-depletable model with no leakage gives
Λ_bg ≈ 0, i.e. today's behavior.

For the numerator there's a check that comes for free. Because k-eff's
denominator is the source weight, which is exactly 1, the reported transport
k-eff is numerically equal to P (all three estimators are estimators of it, so
take the combined value). So production_BOL should reproduce the statepoint
k-eff, and a mismatch means the nu-fission collapse is wrong or the
MicroXS is missing fissile nuclides. Together with §2 covering the
denominator, both terms end up checked.

Assumptions worth naming in the docstring, since they're the ones that degrade
quietly:

• The transport run is the one that produced the fluxes and the MicroXS.
• The fluxes are the ones get_microxs_and_flux returned, in [n-cm/src] per
  source particle, not a renormalized spectrum. Λ_bg is an absolute quantity, so
  unlike the current ratio-only estimator this does not survive rescaling. The
  production_BOL check below catches a violation, and should be an error
  rather than a warning for that reason.
• create_delayed_neutrons matches between that run and the collapse session,
  since Nuclide::nu(E, total) returns prompt-only nu when it's off.
• Non-depletable materials genuinely don't deplete.
• Background absorption per source neutron is constant over the depletion. This
  is the weakest link and it degrades gradually rather than failing loudly.
• All fissile material sits in depletable materials. If you deplete a subset of
  the fuel this is violated, and the natural generalization is a frozen
  background production term symmetric with Λ_bg. Not asking for that here, but
  worth a line in the docstring so nobody is surprised.

If you'd rather not take this on in this PR at all, the alternative is to
restrict the estimate to models where every material is depletable and say so.
I'd rather have the general version.

4. The (n,xn) convention is right, but please show the derivation

I checked this because the appeal to "the definition used elsewhere in OpenMC"
isn't obvious from the Monte Carlo code, where (n,xn) neutrons are physically
created and transported (physics.cpp:1190-1200) and absorption excludes
(n,xn) entirely (endf.cpp:59-70). But writing k-eff in balance form, with the
source weight substituted via 1 = A + L - X:

$$k_\mathrm{eff} = \frac{P}{1} = \frac{P}{A + L - X}$$

That is reduced absorption, which is exactly what OpenMC's k-eff is. Your
treatment is correct. I'd just put those two lines in the docstring so the next
reader can verify it instead of taking it on faith.

I also checked the table programmatically: across all 84 reactions and 351 MTs,
neutrons_out == 0 ⟺ is_disappearance(mt), zero mismatches, including the
easy-to-miss cases like (n,nα) and (n,np) contributing zero. Nice.

5. The volume normalization is missing

This falls out of §1: no definition currently matches the code for more than one
material.

# independent_operator.py:486
rate = atoms * (micro_xs[nuc, rxn] * flux).sum()

self.number[mat, nuc] is total atoms and self.fluxes[i_mat] is
volume-integrated flux in [n-cm/src], so the rate needs the division
_calculate_reaction_rates already applies:

# openmc_operator.py:517-527
volume_b_cm = 1e24 * self.number.get_mat_volume(mat)
...
rates[i] = tally_rates / volume_b_cm

Without it, each material's contribution to both sums carries an extra factor of
V_m, so large-volume regions dominate more than they should. With a single
material it cancels in the ratio, which is why the tests pass. Concretely: a
1 cm³ fuel region plus a 100 cm³ absorber gives the absorber 100 times the
weight it should have, and k comes out far below the true value. Flip the
volumes and the error flips sign.

Since _estimate_k_inf recomputes rates that _calculate_reaction_rates just
produced, one option that fixes this and drops the duplicate Python triple loop
is to derive the loss term from the returned rates array (already per-atom and
volume-normalized), keeping a separate pass only for nu-fission.

Please add a test with two materials of different volumes. All five new tests
use a single material at volume = 1.0, so this class of bug can't fail CI.

6. collapse_nu_fission_rate duplicates Reaction::collapse_rate

src/nuclide.cpp:1075-1162 is essentially a copy of src/reaction.cpp:115-180
with xs_avg replaced by nu(E)*sigma. We already tabulate that product on the
nuclide's union grid:

// src/nuclide.cpp:437-447
xs_[t](i, XS_NU_FISSION) = nu(E, EmissionMode::total) * xs_[t](i, XS_FISSION);

Collapsing that column instead drops the duplicated loop, avoids two virtual
nu() calls per grid segment (order 1e5 segments per nuclide per temperature
for a full-range collapse), and handles evaluations giving partial fission
(MT 19/20/21/38) rather than MT 18. For those, reaction_index_[N_FISSION] is
-1 and the current code silently returns 0.0. Not a regression, since the
existing 'fission' collapse has the same blind spot, but XS_NU_FISSION sums
over all fission reactions so you'd get it for free.

If you'd rather keep the reaction-based path, please factor the shared loop into
a helper taking a per-point weight rather than maintaining two copies.

Thanks for the detailed review. Addressed everything in 3006080:

1. Docstring. Rewrote _estimate_k_inf with the full formula (including 1/V_m), index ranges, the k-inf vs k-eff distinction, and assumptions. Explicitly states the estimate covers depletable materials only and will be higher than the true k-inf when non-depletable absorbers exist.

2. Precondition. After super().__init__() (so self.chain reflects any chain reduction), the operator checks set(self.chain.reactions) <= set(m.reactions) for every MicroXS. Missing channels trigger a warning and disable the estimate. Added test_kinf_disabled_when_reactions_missing.

3. Background absorption. Not implemented in this PR. Documented as a limitation in the docstring (background absorption is assumed constant and ignored).

4. (n,xn) derivation. Added to the docstring: k_eff = P/(A+L-X), reduced absorption A-X matching OpenMC's estimators, done in one pass via the (1-x_j) factor.

5. Volume normalization. Fixed: rate is now divided by volume_b_cm = 1e24 * self.number.get_mat_volume(mat), matching _calculate_reaction_rates. Added test_calculate_kinf_two_materials with three volume pairs (1,4), (2,2), (3,0.5) that fail without the fix.

6. C++ refactor. collapse_nu_fission_rate now collapses the pre-tabulated XS_NU_FISSION column on the nuclide union grid instead of duplicating the integration loop with per-segment nu() calls. No threshold offset needed, and partial fission MTs are included automatically.